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"Universal" Distribution of Inter-Earthquake Times Explained 
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We propose a simple theory for the "universal" scaling law previously reported for the distributions 
of waiting times between earthquakes. It is based on a largely used benchmark model of seismicity, 
which just assumes no difference in the physics of foreshocks, mainshocks and aftershocks. Our 
theoretical calculations provide good fits to the data and show that universality is only approximate. 
We conclude that the distributions of inter-event times do not reveal more information than what 
is already known from the Gutenberg-Richter and the Omori power laws. Our results reinforces the 
view that triggering of earthquakes by other earthquakes is a key physical mechanism to understand 
seismicity. 
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Understanding the space-tinie-magnitude organization 
of earthquakes remains one of the major unsolved prob- 
lem in the physics of the Earth. Earthquakes are char- 
acterized by a wealth of power laws, among them, (i) 
the Gutenberg-Richter distribution ~ 1/ E^^^ (with /? « 
2/3) of earthquake energies E [jj]; (ii) the Omori law 
^ 1/tP (with p Ki 1 for large earthquakes) of the rate of 
aftershocks as a function of time t since a mainshock 2j ; 
(iii) the productivity law ~ E°- (with a < 2/3) giving the 
number of earthquakes triggered by an event of energy 
E |3]; (iv) the power law distribution ^ 1/L^ of fault 
lengths L |j|; (v) the fractal (and evenprobably multi- 
fractal '5'|) structure of fault networks [6^ and of the set 
of earthquake epicenters |3 . The quest to squeeze novel 
information from the observed properties of seismicity 
with ever new ways of looking at the data goes unabated 
in the hope of better understanding the physics of the 
complex solid Earth system. In this vein, from an analy- 
sis of the probability density functions (PDF) of waiting 
times between earthquakes in a hierarchy of spatial do- 
main sizes and magnitudes in Southern California, Bak 
et al. discussed in 2002 a unified scaling law combining 
the Gutenberg-Richter law, the Omori law and the frac- 
tal distribution law in a single framework [8j (see also 
ref. [9'| for a similar earlier study). This global approach 
was later refined and extended by the analysis of many 
different regions of the world by Corral, who proposed the 
existence of a universal scaling law for the PDF H{t) of 
recurrence times (or inter-event times) t between earth- 
quakes in a given region S jlfl lul : 



H{t) ~ A/(At) . 



(1) 



The remarkable finding is that the function f{x), which 
exhibit different power law regimes with cross-overs, is 
found almost the same for many different seismic regions, 
suggesting universality. The specificity of a given region 



seems to be completely captured solely by the average 
rate A of observable events in that region, which fixes the 
only relevant characteristic time 1/A. 

The common interpretation is that the scaling law (Q 
reveals a complex spatio-temporal organization of seis- 
micity, which can be viewed as an intermittent flow of 
energy released within a self-organized (critical?) system 
[l3| , for which concepts and tools from the theory of criti- 
cal phenomena can be applied [13] . Beyond these general 
considerations, there is no theoretical understanding for 
Q. Under very weak and general conditions, Molchan 
proved mathematically that the only possible form for 
fix), if universality holds, is the exponential function 
14', in strong disagreement with observations. Recently, 
from a re-analysis of the seismicity of Southern Califor- 
nia, Molchan and Kronrod [l5J have shown that the uni- 
fied scaling law ^ is incompatible with multifractality 
which seems to offer a better description of the data. 

Here, our goal is to provide a simple theory, which 
clarifies the status of ^, based on a largely studied 
benchmark model of seismicity, called the Epidemic- 
Type Aftershock Sequence (ETAS) model of triggered 
seismicity [laand whose main statistical properties are 
reviewed in [ig. The ETAS model treats all earth- 
quakes on the same footing and there is no distinction 
between foreshocks, mainshocks and aftershocks: each 
earthquake is assumed capable of triggering other earth- 
quakes according to the three basic laws (i-iii) mentioned 
above. The ETAS model assumes that earthquake mag- 
nitudes are statistically independent and drawn from 
the Gutenberg-Richter distribution Q{m). Expressed in 
earthquake magnitudes m oc (2/3) Inio E, the probabil- 
ity Qim) for events magnitudes mi to exceed a given 
value m is Q{m) = io~''(™~"'o), where 6 ~ 1 and mo 
is the smallest magnitude of triggering events. We also 
parametrize the (bare) Omori law for the rate of triggered 



events of first-generation from a given earthquake as 
$(t) = 6'c^/(c+i)^+^ with 9>0. $(t) can be interpreted 
as the PDF of random times of independently occurring 
first-generation aftershocks triggered by some mainshock 
which occurred at the origin of time t = 0. Several au- 
thors have shown that the ETAS model provides a good 
description of many of the regularities of seismicity (see 
for instance Ref. [13 and references therein). 

H{t)/X 




FIG. 1: Fig.l taken from Corral's Ref.^j] plotting the scaled 
(according to (^ probability density functions (PDF) of the 
recurrence times r between successive earthquakes in various 
regions of the world, scaled by their corresponding seismicity 
rates A. Top to bottom: the NEIC worldwide catalog for re- 
gions with L > 180 degrees, 1973-2002; NEIC with L < 90 
degrees, (same period of time); Southern California, 1984- 
2001, 1988-1991, and 1995-1998; Northern California, 1998- 
2002; Japan, 1995-1998, and New Zealand, 1996-2001; (bot- 
tom), Spain, 1993-1997, New Madrid, 1975-2002, and Great 
Britain, 1991-2001. The PDFs have been translated for clar- 
ity. The thin continuous lines are Corral's fits I12II while the 
thick continuous lines are our prediction IIH based on ETAS 
model with the parameters 9 — 0.03, n = 0.9, a — 0.76 and 
p=l. 

Our main result is the theoretical prediction (|ll|l be- 
low, which is used to fit Corral's data in Fig. ^ with re- 
markably good agreement. According to Occam's razor, 
this suggests that the previously mentioned results on 
universal scaling laws of inter-event times do not reveal 
more information that what is already captured by the 
well-known laws (i-iii) of seismicity (Gutenberg-Richter, 
Omori, essentially), together with the assumption that 
all earthquakes are similar (no distinction between fore- 
shocks, mainshocks and aftershocks [l9j), which is the 
key ingredient of the ETAS model. Our theory is able 
to account quantitatively for the empirical power laws 
found by Corral, showing that they result from subtle 
cross-overs rather than being genuine asymptotic scaling 
laws. We also show that universality does not strictly 
hold. 

Our strategy to obtain these results is to first calcu- 



late the PDF of the number of events in finite space-time 
windows 18], using the technology of generating proba- 
bility functions (GPF), which is particularly suitable to 
deal with the ETAS as it is a conditional branching pro- 
cess. We then determine the probability for the absence 
of earthquakes in a given time window from which, us- 
ing the theory of point processes, is determined the PDF 
of inter-event times. Our analysis is based the previous 
calculations of Ref. [ig, which showed that, for large ar- 
eas (L ~ tens of kilometers or more), one may neglect 
the impact of aftershocks triggered by events that oc- 
curred outside the considered space window, while only 
considering the events within the space domain which are 
triggered by sources also within the domain. 

Generating probability functions of the statis- 
tics of event numbers. Consider the statistics of the 
number R(t, r) of events within a time window [i, t + t]. 
It is efficiently described by the method of GPF, defined 
by Qs{z,t) = (z^(*''^)) where the brackets (.) denote a 
statistical average over all possible realizations weighted 
by their corresponding probabilities. We consider a sta- 
tistically stationary process, so that 8s(z, r) does not 
depend on the current time t but only on the window 
duration r. For the ETAS model, statistical stationarity 
is ensured by the two conditions that (i) the branching 
ratio n (or average number of earthquakes/aftershocks 
of first-generation per earthquake) be less than 1 and (ii) 
the average rate to of the Poissonian distribution of spon- 
taneous events henon-zero. The GPF <ds{z,T) can then 
be obtained as [18| 

{z,t) =exp(-iu [l-Q{z,t,T)]dt~iu [1 ~ zQ{z,t)]dt] , 

(2) 

where Q{z, t, t) is the GPF of the number of aftershocks 
triggered inside the window [i, t + t] (i > 0) by a sin- 
gle isolated mainshock which occurred at time and 
<d{z,T) = 0(z,i — 0,t). The first (resp. second) term in 
the exponential in |2Jl describes the contribution of after- 
shocks triggered by spontaneous events occurring before 
(resp. within) the window [i,<-|-r]. 
Ref. [1^] previously showed that 8(2;, i, r) is given by 



e 



e(z,i,T) = G[l-vI/(z,t,r)] , 



(3) 



where G{z) is the GPF of the number of first-generation 
aftershocks triggered by some mainshock, and the auxil- 
iary function ^(z,i,T) satisfies to 

*(z, i, t) = [1 - e(z, t, T)]®^{t)+[l ~ zQ{z, r)]®$(i+T) . 

(4) 
The symbol ® denotes the convolution operator. Inte- 
grating Q with respect to t yields J^ 'i{z,t,T) dt ~ 
J^[l - Q{z,t,T)]dt+ [1 - zQ{z,t)] iS>a{T), so that ex- 
pression (|2Jl becomes 

/■ 00 

es(2:,T)=exp -iu *(2:,t,r)dt-tj[l - 2:6(2, r)] (g)&(r) 

L Jo 

(5) 

where b{t) = /J $(t') dt' and a{t) = 1 - b{t) ^ j^^- 



Probability of abscence of events. For our pur- 
pose, the probability PsIt) that there are no earthquakes 
in a given time window of duration r provides an intuitive 
and powerful approach. It is given by 

/■oo 

Po(r) = 65(2: = 0,r) =exp -uj *(i, r) dt - a;r + Ci;A(r) 

L Jo 

(6) 

where *(i,T) = *(z = 0,i, r) and A{t) = J^ a{t)dt ~ 

TTe (T/c)^"^ for T>c). ^^ 

To make progress in solving (|3I4I5|I . let us expand G(z) 
in powers of z: 



G{z) ^ I - n + nz + (3{1 - z^ + . 



(7) 



where 7 = b/a (where a = (3/2)a < 1 is the productivity 
exponent when using magnitudes) and /3 = nr(2 — 7)(7 — 
i)'^^^/Y'- While we can calculate the looked-for distri- 
bution of recurrence times using the shown expansion up 
to order (1 — z)^ , it turns out that truncating ^ at the 
linear order is sufficient to explain quantitatively Corral's 
results, as we show below. Using G{z) = 1 — n + nz has 
the physical meaning that each earthquake is supposed to 
generate at most one first-generation event (which does 
not prevent it from having many aftershocks when sum- 
ming over all generations). Indeed, interpreted in proba- 
bilistic term, G(z) = 1 — ji + nz says that any earthquake 
has the probability 1 — 71 to give no offspring and the 
probability n to give one aftershock (of first-generation). 
This linear approximation is bound to fail for small re- 
currence times associated with the large productivity of 
big earthquakes and, indeed, we observe some deviations 
for the shorter recurrence times below several minutes 
as discussed below. The linear approximation is not in- 
tended to describe the statistics of very small recurrence 
times within clusters of events triggered by large main- 
shocks, but is approprioate for "quiet" periods of seismic 
activity. The linear approximation is expected and ac- 
tually seen to work remarkably well for large recurrence 
times of hours, days, weeks... 

The linear approximation bypasses much of the com- 
plexity of the nonlinear integral equations (|3I4|) to obtain 
J^ ^(i, T)dt = j^- Expression IjHJl becomes (for t ^ c) 



Ps (r) = exp 



-(1 



nap 
T~6 



(8) 



where 



x = \t, a=(Aoc)^ p^A/Ao = g(TO) ( — j .(9) 

The average seismicity rate A is given by A = j^ , which 
renormalizes the average rate u; of spontaneous sources 
by taking into account earthquakes of all generations trig- 
gered by a given source: A = w + ncj + n^uj + .... Due to 
the assumed statistical independence between event mag- 
nitudes, the proportion between spontaneous observable 



events and their observable aftershocks does not depend 
on the magnitude threshold and the above expression for 
the average seismic rate holds also for observable events 
at different magnitude thresholds of completeness. Fi- 
nally, Ao is the average seismic rate within a spatial do- 
main So of reference with linear size Lq, and p takes into 
account the dependence on the magnitude threshold m 
for observable events and on the scale L of the spatial 
domain S used in the analysis. The first term (1 — n)x 
in the exponential of JHJ describes the exponential de- 
creasing probability of having no events as r increases 
due to the spontaneous occurrence of sources. The other 
term proportional to x^~^ takes into acount the influ- 
ence through Omori's law of earthquakes that happened 
before the time window. 

Statistics of recurrence times. Consider a sequence 
of times {U} of observable earthquakes, occurring inside 
a given seismic area iS. The inter-event times are by def- 
inition Ti = ti — U^i. The whole justification for the cal- 
culation of Ps (r) lies in the well-known fact in the theory 
of point processes |2flj that the PDF H{t) of recurrence 
times Ti is given by the exact relation 



H{t) 



1 d^Psjr) 

A dT^ 



(10) 



Substituting IjHJl in this expression yields our main theo- 
retical prediction for the PDF of recurrence times, which 
is found to take the form Q with 



f{x) = (anOp^x^'^^^ +[l-n + nap'^x^'^]^'] 
exp —(1 — n)x X 



(11) 



While our theoretical derivati on j ustifies the scaling re- 
lation (0) observed empirically [llJ, lll| , the scaling func- 
tion f{x) given by ((TT1) is predicted to depend on the 
criticality parameter n, the Omori law exponent 9, the 
detection threshold magnitude m and the size L of the 
spatial domain S under study. While 6 might perhaps 
be argued to be universal, this is less clear for n which 
could depend on the regional tectonic context. The sit- 
uation seems much worse for universality with respect 
to the two other parameters m and L which are catalog 
specific. It thus seems that our prediction can not agree 
with the finding that f{x) is reasonably universal over 
different reg ions of the world as well as for worldwide 
catalogs lldllll]. 

It turns out that the dependence on the idiosynchratic 
catalog-dependent parameters m and L is basically irrel- 
evant as long as 9 is small and n in the range 0.7 — 1 
previously found to be consistent with several other sta- 
tistical properties of seismicity [l^, |2l|. Note that the 
condition that 9 be small is fully compatible with many 
empirical studies in the literature for the Omori law re- 
porting an observable (renormalized) Omori law decay 



- l/t°-9-i corresponding to = 0-0.1 [13. Fig. H 
shows the changes of f{x) when varying the magnitude 
threshold from to 3. These changes of f{x) seem to 
be within the inherent statistical uncertainties observed 
in empirical studies |10l lll| . The technical origin of the 
robustness lies in the fact that, for 9 = 0.03 say, chang- 
ing m — mo from to 6 amounts to changing p from 1 
(to = TTio) to p ~ 10^^ {m = Too + 6) which changes p^ 
from 1 to only p^ ~ 0.66. We conclude that our theory 
provides an explanation for both the scaling ansatz (^ 
and its apparent universal scaling function. 




FIG. 2: Scaling function f{x) defined in l|ll|l for n — 0.8, 
9 = 0.03, a = 0.76 and for two values of m — mo = 0, 3 
corresponding to a thousand-fold variation of p = 1; 10""^. In 
these synthetic examples, we assume that the spatial domain 
So corresponds to an average seismicity rate Ao — 1 per day, 
that the characteristic time scale of the Omori law is c ~ 10 
sec, so that Ac — 10~^. Then, for 6 ~ 0.03, we have a ~ 
0.76. The obtained function f{x) is compared with Corral's 
empirical fitting function fc{x) defined in 11211 with g = 0.7, 
S = 1.05, d = 1.7 and C = 0.78. The dashed lines are the 
power laws ~ x'"'^ and x~^ . 

We can squeeze more out of Ijllf) to rationalize the 
empirical power laws reported by Corral. In particular. 
Corral proposed the following empirical form for f{x) 
which, in our notations, reads 



fc{x) 



C6 



drh/s) 



Q 



9-1 



,-Wd)' 



(12) 



where g = 0.67±0.05, 6 = 1.05±0.05, d = 1.64±0.15, and 
C ensures normalization [l0|,lll|- Fig|21shows indeed that 
expression (I12II with Corral's reported parameter values 
for g, 5 and d fits 1)11(1 remarkably well quantitatively. 
While the intermediate asymptotics f{x) ~ x^~'^ ~ x~'^^^ 
proposed by Corral is absent from our theoretical expres- 
sion Hll() . it can actually be seen as a long cross-over 
between the power and exponential factors in ((ll|l . as 
shown by one of the dashed lines in Fig [3 

Interestingly, expressions Hll|) and (|12|) depart from 
each other for x < 0.01. Our theoretical distribution 



f{x) has the power law asymptotic f{x) ^ x~^, which is 
a direct consequence of Omori's law described explicitly 
by the first power law factor in front of the exponential 
in Hill) . It is absent from expression (|12|l . However, its 
presence is clear in real data as shown in Fig. ^ extracted 
from |ll| on which we have superimposed our theoreti- 
cal prediction (|ll|l . Note that expression Hll() exhibits 
a slight departure from the data for small x's (defined 
in ©), which can be attributed to the linearization of 
(|7|). which amounts to neglecting the renormalization of 
the Omori law by the cascade of triggered aftershocks 
[ITJ. Taking into account this renormalization effect by 
the higher-order terms in the expansion (Q improves the 
fit to the data shown in Fig. ^ Our detailed study shows 
that comparing (|ll|l with data provides constraints on 
the parameter n: the data definitely excludes small val- 
ues of n and seems best compatible with n = 0.7 — 1, in 
agreement with previous constraints [l^l suggesting that 
earthquake triggering is a dominant process. 
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